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Introduction 

There  is  currently  much  interest  in  the  unsteady  aerodynamics  of  flows  relevant  to  rotorcraft 
(manned  and  unmanned)  air  vehicles.  The  flow  fields  in  these  applications  exhibit  unsteady 
separation  followed  by  the  formation  of  dynamic-stall-like  vortices  whose  evolution  and 
interaction  with  flexible  aerodynamic  surfaces  have  a  significant  impact  on  flight  stability  and 
performance.  Analysis  of  these  flows  is  complicated  further  by  their  mixed  laminar  transitional 
turbulent  character  at  moderate  Reynolds  numbers,  as  well  as  by  the  broad  range  of  possible 
parameters,  kinematics,  and  configurations.  To  facilitate  progress  in  the  understanding  and 
prediction  of  the  relevant  fluid  dynamic  mechanisms,  it  is  natural  to  consider  methods  that 
provide  simplification  of  the  flow  phenomena  by  separating  them  into  individual  modes.  The 
technique  of  Proper  Orthogonal  Decomposition  (POD)  is  a  popular  way  of  accomplishing  this 
task.  However,  while  POD  is  capable  of  extracting  the  most  energetic  parts  of  the  flow  field,  it 
has  been  shown  to  lack  ability  of  highlighting  subtle,  oscillatory  phenomena  that  are  nevertheless 
implicated  in  important  physical  processes  such  as  the  shear  layer  separation. 

Unsteady  flows  over  plunging  and  pitching  airfoils  with  large  excursions  in  effective  angle  of 
attack  exhibit  the  phenomenon  termed  dynamic  stall,  a  process  characterized  by  unsteady 
separation  and  by  the  formation  of  large-scale  leading-edge  and  trailing-edge  vortices,  which 
exert  difficult-to-predict  variations  in  aerodynamic  loads.  Comprehensive  reviews  of  this 
phenomenon,  first  discovered  and  studied  extensively  in  the  context  of  helicopter  rotor  blades, 
has  been  given  in  [1,3,8]. 

The  initial  research  has  suggested  existence  of  high-frequency  effects  in  leading  edge  vortex 
shedding  events,  not  present  in  the  static  stall  case.  The  purpose  of  the  project  was  to  advance 
this  approach  using  Aimdyn's  development  of  Koopman  operator  mode  decomposition 
techniques. 
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Computation  of  Koopman  Modes 

Koopman  mode  decomposition  is  based  on  the  surprising  fact,  discovered  in  [5],  that  normal 
modes  of  linear  oscillations  have  its  natural  analogue  -  Koopman  modes  -  in  the  context  of 
nonlinear  dynamics.  To  pursue  this  analogy,  one  must  change  the  representation  of  the  system 
from  the  state-space  representation  to  the  dynamics  governed  by  the  linear  Koopman  operator 
([4])  on  an  infinite-dimensional  space  of  observables.  Contrary  to  the  proper  orthogonal 
decomposition,  the  dynamic  mode  decomposition  contains  not  only  information  about  coherent 
structures,  but  also  about  their  temporal  evolution. 

Based  on  snapshots  of  the  flow,  we  can  approximate  the  Koopman  Modes  using  an  Amoldi-like 
algorithm  sometimes  called  dynamic  mode  decomposition  (DMD)  ([9,10])  which  computes 
eigenvalues  based  on  the  so-called  companion  matrix. 


Given  a  sequence  of  equispaced  in  time  snapshots  from  numerical  simulations  or  physical 
experiments,  with  At  being  the  time  interval  between  snapshots,  a  data  matrix  is  formed  with 
columns  that  represent  the  individual  data  samples  Uj  G  R^,j  =  0, ... ,  m,  with  j  representing  time 
jAt.  The  companion  matrix  is  then  defined  as: 


C  = 


/  0  0  •••  0 

1  0  0 

0  1  0 


CO  \ 

Cl 

C2 


y  0  0  •  •  •  1  Cjn—i  j 


where  Cf,  i  =  0, ... ,  m  —  1  are  such  that: 

m-l 

Ujn  ~  CjU|  -j-  r 
J=0 


and  r  is  the  residual  vector. 


The  spectrum  of  the  Koopman  operator  restricted  to  the  subspace  spanned  by  Uj  is  equal  to  the 
spectrum  of  the  infinite-dimensional  companion  matrix  and  the  associated  Koopman  modes  are 
given  by  Ka  (provided  that  a  does  not  belong  to  the  null  space  of  K),  where 
K  =  {uq,u^,  ...  ,Ufn-i\  is  the  column  matrix  (vector-valued)  of  observables  snapshots  at  times 
0,  At, ...,  (m  —  l)At  and  a  is  an  eigenvector  of  the  shift  operator  restricted  to  Krylov  subspace 
spanned  by  Uj  which  the  Companion  matrix  is  an  approximation  of  The  approximate  Koopman 
eigenvalues  and  eigenvectors  obtained  by  the  Amoldi's  algorithm  are  called  Ritz  eigenvalues  and 
eigenvectors. 

The  standard  Amoldi-type  algorithm  to  calculate  the  Ritz  eigenvalues  Ay  and  eigenfunctions  Vy  is 
as  follows: 

1.  Define  K  =  [iio,Ui,  ...,u^_i]. 

2.  Find  constants  Cy  such  that: 
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m— 1 

»•  =  «r.!  -  ^  CjUj  =  Urn  “  A'c,  r±span{uQ, Um-i}- 
j=o 

This  can  be  done  by  defining  c  =  where  K'^  is  the  pseudo  inverse  of  K. 


3.  Define  the  companion  matrix  C  as  above  and  find  its  eigenvalues  and  eigenvectors: 

C  =  T~^AT,A  =  diag{Xi,-  ■  ■  ,  A„,), 

where  eigenvectors  are  columns  of  7“^.  Note  that  the  Vandermonde  matrix  f 

m—1 


/  1  Ai  A? 

1  A2  A| 


V  1 


A?-* 


/ 

diagonalizes  the  companion  matrix  C  ,  as  long  as  the  eigenvalues  are  distinct. 


4. 


Define  Vj  to  be  the  columns  of  K  =  Kf 


Then,  the  Amoldi-type  Koopman  Mode  Decomposition  gives: 

m 

VA-  =  [0, 1,  •  •  •  ,  m],  Uk  =  >>jV{:,j) 

j=l 


In  order  for  the  result  of  an  Amoldi-type  method  to  give  a  good  approximation  of  the  Koopman 
modes,  the  companion  matrix  C  should  be  a  good  matrix  representation  of  the  projection  of  the 
Koopman  operator  on  a  Krylov  subspace.  The  algorithm  described  above  uses  a  least-square 
approximation.  However,  the  known  properties  of  the  Koopman  operator  are  ignored  in  the 
above  algorithm: 

■V  Constant  functions  are  eigenfunctions  of  the  Koopman  operator  at  eigenvalue  1 .  Hence  the 
Companion  matrix  as  an  approximation  of  the  Koopman  operator  should  have  an  eigenvalue 
at  1.  This  condition  can  be  translated  as:  Cj  =  1. 

V  When  the  signal  is  periodic  on  the  attractor  such  that  Uq  =  u^,  the  DMD  computation  should 
reduce  to  Discrete  Fourier  Transform  [2].  In  that  case,  the  companion  matrix  coefficient 
should  be  Cq  =  l,Cj  =  0,  V;  >  0  and  the  algorithm  computing  the  Koopman  Modes  should 
be  able  to  approximate  those  values  of  Cj  automatically. 


A  new  Amoldi-type  algorithm  has  been  developed  at  Aimdyn.  It  is  similar  to  the  standard  one 
but  includes  the  known  properties  described  above  by  replacing  step  2  of  the  algorithm  where  the 
constants  Cj  are  computed  by: 

c  =  K+Um, 

where 


and  is  the  pseudo  inverse  of  K. 


U2 

1 


Um-{  \  f  Um-  m  \ 

1  =  [  0  j 
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Analysis  of  a  2D  simulation  data  for  1  period  of  airfoil  oscillation 

By  using  the  described  above  techniques  a  simulation  data  set  from  Dr.  Miguel  Visbal  at  the 
AFRL  was  analyzed.  The  files  were  360  snapshots  of  2D  grid  and  flow  data  for  1  period  of 
airfoil  oscillation.  The  flow  was  spanwise-averaged  so  that  small  scales  structures  are  reduced. 
The  obtained  data  was  processed  for  plotting  (see  Figure  1  that  shows  u  horizontal  velocity,  v 
vertical  velocity,  w  out  of  plane  velocity,  p  pressure  and  rho  density  at  phase  %I2). 


Figure  1.  Dynamic  Stall  of  a  Plunging  Airfoil  2D  simulation  data  sets  at  phase  nil  of  the 
plunging  motion,  u  horizontal  velocity,  v  vertical  velocity,  w  out  of  plane  velocity,  p  pressure 

and  rho  density. 

Due  to  the  large  size  of  the  data  set,  using  Matlab  for  data  processing  reaches  its  limit.  AIMdyn’s 
team  implemented  the  Koopman  Mode  Decomposition  algorithm  in  C++  to  be  able  to  process 
very  large  sets  of  data. 

A  new,  Amoldi-type  method  developed  by  Aimdyn  was  used  for  Koopman  Mode 
Decomposition.  Figure  2  shows  Fourier  and  KMD  spectrum  for  u-velocity. 


Figure  2.  Fourier  and  KMD  spectrum  for  u-velocity.  (left)  Frequencies  from  0  to  0.5.  (right) 
Close-up  for  frequencies  from  0  to  0.04.  The  amplitude  is  calculated  as  the  norm  of  all 
amplitudes  over  all  spatial  points.  FFT  is  in  black;  unstable  KM  spectrum  is  red;  stable  KM 

spectrum  is  green. 
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The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  3. 


Eigenvalues  on  the  Unit  Circle 


Figure  3.  Koopman  Mode  Eigenvalues. 


Fourier  and  KMD  spectrum  for  several  regions  of  interest  for  u  velocity  (see  Figure  4)  are  shown 
in  Figure  5. 


1.5 

1 

0.5 

0 


Figure  4.  Regions  of  interest  for  u  velocity. 


u-velocity 


R1  R2  R3  R4 


Figure  5.  Fourier  and  KMD  spectrum  for  four  regions  of  interest  for  u-velocity.  FFT  is  in  black; 
unstable  KM  spectrum  is  red;  stable  KM  spectrum  is  green. 

The  frequency  spectrum  and  the  magnitude  of  modes  corresponding  to  frequency  0.027  Hz  and 
frequency  0.01 1  Hz  for  the  u-velocity  probe  with  coordinates  (1.04;  0.50)  are  shown  in  Figure  6. 
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Figure  6.  Fourier  and  KMD  spectrum  for  a  point  of  u-velocity  with  (1.04;  0.50)  coordinates.  FFT 
is  in  black;  KM  spectrum  is  blue.  Magnitude  for  the  Koopman  mode  corresponding  to  frequency 

0.027  Hz  and  frequency  0.01 1  Hz. 

Using  Koopman  mode  decomposition  technique  processing  of  a  2D  data  set  the  team  members 
studied  details  of  the  physics  of  dynamic  stall. 


Similarity  of  the  Koopman  mode  decomposition  of  the  dynamic  stall  and  the  cylinder  wake 
with  oscillating  Re  number  forcing 

We  developed  a  model  system  that  is  capable  of  revealing  physics  of  the  d5mamic  stall  in  a 
simpler  setting.  The  model  system  is  a  cylinder  in  an  incoming  oscillatory  flow.  It  provided  us 
with  some  remarkable  insight.  In  the  following  we  describe  the  simulation  and  some  results. 

The  cylinder  wake  simulation  was  provided  by  Dr.  Bryan  Glaz  at  the  U.S.  Army  Research 
Laboratory.  The  cylinder  was  0.002m  in  diameter.  Flow  snapshots  were  outputted  starting  at  the 
25th  time  step,  up  to  the  60,000th  time  step  in  increments  of  25  time  steps  (each  time  step  is  le-4 
s).  The  total  time  length  corresponds  to  12  periods  of  the  driving  frequency,  which  is  2Hz.  The 
number  of  nodes  in  the  grid  is  29954. 

From  0.0s  to  0.5s  (i.e.  5000  iterations)  Reynolds  number  (Re)  was  58.3  that  corresponds  to 
incoming  velocity  of  0.5  m/s.  The  critical  Re  for  a  cylinder  is  about  Re~40  when  it  starts 
exhibiting  the  von  Karman  wake  instability.  So  the  interval  from  0  to  0.5  s  at  Re=58.3 
corresponds  to  the  Landau  equation  when  the  external  input  is  greater  than  the  bifurcation  value. 
Then,  starting  at  0.5  seconds,  the  Re  is  being  oscillated  by  oscillating  incoming  velocity.  The 
oscillating  Re  corresponded  to  Re  =  58.3  +  35*sin(2pi*omega*t),  where  omega  corresponded  to 
2  Hz.  The  reason  35  was  selected  for  the  oscillating  Re  amplitude  is  because  a  large  enough 
amplitude  was  needed  such  that  the  Re  would  oscillate  above  and  below  the  critical  value 
predicted  by  the  regular  Landau  equation.  The  reason  2  Hz  was  selected  is  that  it's  an  order  of 
magnitude  slower  than  the  von  Karman  vortex  frequency,  which  is  about  39  Hz  for  Re  =  58.3.  It 
was  planned  so  that  the  driving  frequency  to  be  about  an  order  of  magnitude  lower  than  the  flow 
instability  frequency  because  the  dynamic  stall  driving  frequency  is  about  1-2  orders  of 
magnitude  slower  than  the  shear  layer,  wake  dynamics  for  the  airfoil.  So  the  goal  was  to  setup  a 
cylinder  problem  where  there  is  a  significant  separation  in  time  scales. 
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The  Koopman  spectram  plots  are  shown  in  Figure  7.  It  is  remarkable  how  the  spectrum  shown  at 
the  bottom  becomes  broad  upon  introduction  of  the  low-frequency  oscillation,  in  contrast  to  the 
clean  peak  at  40  for  the  von  Karman  vortex  flow  in  the  top  of  that  figure.  An  analytical 
understanding  of  that  is  lacking,  despite  a  number  of  dynamical  system  bifurcation  studies  (e.g. 
[7])  and  we  plan  to  investigate  it  in  the  future  work. 
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Figure  7.  Spectrum  of  von  Karman  vortex  shedding  flow  around  a  cylinder  (top)  and  the  flow 
around  the  same  cylinder  when  the  incoming  velocity  is  oscillated  (bottom).  Note  the  broadening 
of  the  spectrum,  similar  to  the  pitching  wing  spectrum,  and  the  high  frequency  component 

around  120. 


The  interesting  spectrum-broadening  effect  for  flow  around  cylinder  in  incoming  flow  of 
oscillating  magnitude  leads  us  to  believe  that  further  investigation  of  this  model  might  shed 
physical  insight  on  the  wing  pitching  and  plunging  model,  that  also  has  broad  spectrum. 

In  [2],  the  observation  is  made  that  subtracting  the  mean  of  the  sequence  of  snapshots  leads  to 
the  result  of  all  possible  eigenvalues  being  on  the  unit  circle,  the  companion  matrix  analysis 
reducing  essentially  to  the  Discrete  Fourier  Transform,  however  in  [6],  it  is  shown  that  this  is 
true  for  finite  m  only  if  the  observable  snapshots  are  periodic  with  Uq  = 

For  the  described  above  cylinder  wake  simulation  the  Koopman  Mode  Decomposition  (by  using 

Amoldi-type  algorithm)  was  performed  for: 

case  A:  subtracting  the  mean  of  the  sequence  of  snapshots; 

case  B:  without  subtracting  the  mean  of  the  sequence  of  snapshots. 

Case  B  was  studied  for  the  two  following  sub-cases: 

case  Bl:  not-normalizing  each  column  of  the  inverse  to  the  Vandermonde  matrix 
case  B2:  normalizing  each  column  of  the  inverse  to  the  Vandermonde  matrix  f~^. 
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Note  that  in  the  case  of  subtracting  the  mean,  the  normalized  and  not-normalized  cases  are  very 
similar. 


Figure  8  shows  Fourier  and  KMD  spectrum  for  u-velocity  for  the  whole  simulation  time  for 
cases  A,  B1  and  B2. 
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Figure  8.  Fourier  and  KMD  spectrum  for  u-velocity  for  the  whole  simulation  time.  FFT  is  in 

black;  KM  spectrum  is  red. 
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The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  9. 
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Figure  9.  Koopman  Mode  Eigenvalues  for  u-veloeity  for  the  whole  simulation  time. 
Figure  10  shows  the  DMD  speetrum  in  frequeney  w,  exponential  rate  mu  plane. 

ease  A  easeBl  easeB2 


Figure  10.  DMD  speetrum  in  frequeney  w,  exponential  rate  mu  plane  for  u- velocity  for  the 

whole  simulation  time. 


Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  modes  corresponding  to 
frequency  2  and  63.83  Hz  for  case  A  are  shown  in  Figure  1 1 . 
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Figure  11.  Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  mode  for  u-veloeity 
for  the  whole  simulation  time  for  case  A  corresponding  to  frequency  2  Hz  (left)  and  63.83  Hz 

(right). 


The  recomposition  of  signal  using  20  pairs  of  the  highest  magnitude  modes  sorted  by  abs  of 
norm  of  Vj  has  been  performed.  Figure  12  shows  the  reconstructed  signal  (red)  vs  the  original 
data  (blue)  for  the  random  location  and  the  dependence  of  the  mean  of  the  recomposition  error 
(L2  norm  of  the  difference  between  signal  and  recomposed  signal  divided  by  the  L2  norm  of  the 
signal)  on  the  number  of  modes  used. 

case  A  caseBl  caseB2 
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(Bottom)  the  dependence  of  the  mean  of  the  recomposition  error  (L2  norm  of  the  difference 
between  signal  and  recomposed  signal  divided  by  the  L2  norm  of  the  signal)  on  the  number  of 
modes  used  for  u  velocity  for  the  whole  simulation  time. 


10 


Contract  Number:  W911NF-13-C-0081 
Proposal  Number:  64018EG 


Figure  13  shows  the  histogram  of  the  recomposition  error  over  all  locations  if  using  only  4 
modes  (frequencies  2  Hz,  -2  Hz,  63.83  Hz  and  -63.83  Hz)  and  the  recomposition  of  signal  that 


i 
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Figure  13.  (Left)  The  histogram  of  the  recomposition  error  over  all  locations  if  using  only  4 
modes  (frequency  2,  -2,  63.83  and  -63.83),  (Right)  the  recomposition  of  signal  (red  - 
reconstructed  signal,  blue  -  data)  that  gives  the  median  error  for  u  velocity  for  the  whole 

simulation  time  for  case  A. 


Figure  14  shows  Fourier  and  KMD  spectrum  for  u- velocity  for  the  first  period  of  the  simulation 
time  after  removing  first  500  iterations  for  cases  A,  B1  and  B2. 
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Figure  14.  Fourier  and  KMD  spectrum  for  u- velocity  for  the  first  period  of  the  simulation  time. 

FFT  is  in  black;  KM  spectrum  is  red. 
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The  obtained  Koopman  mode  eigenvalues  are  shown  in  Figure  15. 
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Figure  15.  Koopman  Mode  Eigenvalues  for  u- velocity  for  the  first  period  of  the  simulation  time. 
Figure  16  shows  the  DMD  spectrum  in  frequency  w,  exponential  rate  mu  plane. 

case  A  caseBl  case  B2 


Figure  16.  DMD  spectrum  in  frequency  w,  exponential  rate  mu  plane  for  u-velocity  for  the  first 

period  of  the  simulation  time. 

Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  modes  for  case  A  corresponding 
to  frequency  2  Hz  and  38  Hz  are  shown  in  Figure  17. 
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Figure  17.  Magnitude,  phase,  real  part  and  imaginary  part  for  the  Koopman  mode  for  u-veloeity 
for  the  first  period  of  the  simulation  time  for  case  A  corresponding  to  frequency  2  Hz  (left)  and 

38  Hz  (right). 

The  recomposition  of  signal  using  20  pairs  of  the  highest  magnitude  modes  sorted  by  abs  of 
norm  of  Vj  has  been  performed.  Figure  18  shows  the  reconstructed  signal  (red)  vs  the  original 
data  (blue)  for  the  random  location. 


case  A  caseBl  caseB2 


Figure  18.  The  reconstructed  signal  (red)  vs  the  original  data  (blue)  for  a  random  location  for  u 
velocity  for  the  first  period  of  the  simulation  time. 

The  Koopman  Mode  Decomposition  of  the  model  system  that  we  have  pursued  shows 
remarkable  ability  to  reconstruct  signal  coming  from  such  a  broad-spectrum  flow. 
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Conclusions 

The  initial  analysis  of  dynamic  stall  using  the  methods  of  Koopman  operator  theory  was 
performed.  This  research  guided  by  dynamical  systems  methods  and  efficient  numerical  methods 
for  Koopman  Mode  Decomposition  lead  to  some  interesting  outcomes.  Among  those  are  insight 
into  the  physical  effects  of  the  oscillatory  nature  of  the  incoming  velocity  that  the  wing 
experiences  due  to  pitch  or  plunging.  For  example,  the  examination  of  a  model  system  -  a 
cylinder  in  an  incoming  oscillatory  flow  -  showed  some  of  the  physical  effects  (e.g.  broadening 
of  the  spectrum)  observed  in  the  pitching  airfoil  case  are  present  in  the  model  system  and  shed 
light  on  the  dynamics  of  aspects  of  pitching  airfoil  dynamics  as  a  nonlinear  interaction  of  the 
forcing  by  the  oscillating  incoming  flow  frequency  and  natural  frequency  of  vortex  shedding 
dynamics. 
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